/*
 * EquilibriumSolver.h
 *
 *  Created on: 26 Jul 2011
 *      Author: Allan
 */

#ifndef EQUILIBRIUMSOLVER_H_
#define EQUILIBRIUMSOLVER_H_

// C++ includes
#include <vector>
#include <functional>
using namespace std;

// Eigen includes
#include <Eigen/Core>
using namespace Eigen;

// GeoReact includes
#include "ReactionSystem.h"
#include "EquilibriumConditions.h"
#include "EquilibriumProblem.h"
#include "Multiphase.h"

#include "Numerics/MultirootSolver.h"

class EquilibriumSolver 
{
public:
	struct Diagnostic
	{ unsigned iterations; bool converged; };
	
public:
	EquilibriumSolver();
	
	EquilibriumSolver(const Multiphase& multiphase, const ReactionSystem& reactions);
	
	void Initialiaze(const Multiphase& multiphase, const ReactionSystem& reactions);
	
	void SetEquilibriumConditions(const vector<EquilibriumCondition>& eConditions);
	
	void SetSolverOptions(double tolerance, unsigned maxiterations);
	
	const Diagnostic Solve(double T, double P, VectorXd& n);
	
private:
	/// The multiphase instance
	Multiphase multiphase;
	
	/// The system of equilibrium-controlled reactions in the canonical representation
	ReactionSystem eReactions;
	
	/// The equilibrium problem object which overloads the methods Function and Jacobian from MultirootProblem 
	EquilibriumProblem problem;
	
	/// The multiroot solver object used to solve the non-linear system of equations
	MultirootSolver solver;
	
	/// The tolerance for the iterative procedure
	double tolerance;
	
	/// The maximum number of iterations for the iterative procedure
	unsigned maxiterations;
	
private:
	void UnderRelax(VectorXd& n, const VectorXd& nOld) const;
};

#endif /* EQUILIBRIUMSOLVER_H_ */
